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Abstract 

The dynamic critical behavior of isotropic Heisenberg ferromagnets with a planar free surface is 
investigated by means of field-theoretic renormalization group techniques and high-precision com- 
puter simulations. An appropriate semi- infinite extension of the stochastic model J is constructed. 
The relevant boundary terms of the action of the associated dynamic field theory are identified, 
the implied boundary conditions are derived, and the renormalization of the model in d < 6 bulk 
dimensions is clarified. Two distinct renormalization schemes are utilized. The first is a massless 
one based on minimal subtraction of dimensional poles and the dimensionality expansion about 
d = 6. To overcome its problems in going below d = 4 dimensions, a massive one for fixed dimen- 
sions d < 4 is constructed. The resulting renormalization group (or Callan Symanzik) equations 
are exploited to obtain the scaling forms of surface quantities like the dynamic structure factor. 
In conjunction with boundary operator expansions scaling relations follow that relate the critical 
indices of the dynamic and static infrared singularities of surface quantities to familiar static bulk 
and surface exponents. To test the predicted scaling forms and scaling-law expressions for the crit- 
ical exponents involved, accurate computer-simulation data are presented for the dynamic surface 
structure factor. These are in conformity with our predictions. 

PACS numbers: 75.10.Hk, 68.35.Rh, 64.60.Ht, 05.70. Jk 
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I. INTRODUCTION 



A cornerstone of the modern theory of critical phenomena is the arrangement of mi- 
croscopically different systems in universality classes of equivalent critical behavior.0! Few 
basic properties, such as the spatial dimension d, the order-parameter symmetry, and gross 
features of the interactions determine to which universality class for static bulk critical be- 
havior a particular system belongs. These universality classes can be represented by simple 
continuum models like the 4 model, which are minimal in the sense that dropping any of 
the Hamiltonian's terms implies a change of the universality class. An important alternative 
way of representing the universality classes is through standard lattice (spin) models such 
as the Ising model, which lend themselves best for precise Monte Carlo simulations. 

A similar classification scheme exists for dynamic bulk critical behavior.! The associ- 
ated universality classes — called dynamic bulk universality classes henceforth — additionally 
depend on basic properties of the dynamics such as conservation laws, and since distinct dy- 
namics may have the same equilibrium distribution, each static universality class generally 
splits up into several dynamic ones. The latter are represented by stochastic models called 
A, B, ... , J.i 

Research over the past 25 years has revealed the existence of analogous universality classes 
for static surface critical behavior of semi-infinite systems at bulk critical points^ To which 
static surface universality class a given system belongs is decided by its static bulk univer- 
sality class and additional relevant surface properties. Hence each static surface universality 
class as well as each dynamic bulk universality class usually splits up into separate dynamic 
surface universality classes. Furthermore, systems belonging to the same static surface uni- 
versality class and the same dynamic bulk universality class may be representative of distinct 
dynamic surface universality classes as local changes of the dynamics at the surface can be 
relevant .H!! 

Unfortunately, the number of detailed theoretical investigations of dynamic surface crit- 
ical behavior performed until now is rather limited.0!'!00lll Furthermore, they focused 
more or less exclusively on models with purely relaxational dynamics. On the experimental 
side, the situation is worse: stringent experimental checks of the theoretical predictions for 
dynamic surface critical behavior, though urgently needed, are still lacking. One obvious 
reason for this is the difficulty of such experiments. The impressive progress made during 
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the past two decades in the perfection of surface-sensitive scattering techniques has so far 
led only to accurate experimental investigations of static surface critical behavior.0'@Ellll 
Demonstrating that similarly conclusive data can also be obtained for dynamic surface crit- 
ical behavior remains a major experimental challenge, albeit such experiments are expected 
to become feasible in the near future. According to the recent TESLA design report ,0 the 
X-ray free electron laser offers a great potential for such experiments. 

Theoretical progress can play an important role in stimulating such experiments. We 
believe that theoretical advances in two directions are essential for achieving this goal. On the 
one hand, models representing other bulk dynamic universality classes must be considered, 
generalized to systems with boundaries, and carefully investigated to find out what kinds 
of dynamic surface critical behavior can occur, i.e., which dynamic surface universality 
classes exist. On the other hand, detailed theoretical predictions should be worked out for 
experimentally accessible quantities like structure functions etc. 

Pursuing these goals, we will investigate the dynamic surface critical behavior of isotropic 
Heisenberg ferromagnets in this paper. Well-known characteristic features of the dynamics of 
such magnets are the presence of nondissipative (mode- coupling) terms and the conservation 
of the order parameter. We shall employ two different lines of approaches: (i) analytic work 
based on the field-theoretic renormalization group (RG) and (ii) computer-simulation studies 
of the dynamic surface structure function. A brief account of parts of our work has been 
given elsewhere. @ 

In our RG work we utilize an appropriate semi-infinite extension of the usual stochastic 
bulk model J,!Sll^H0iil which represents the dynamic bulk universality class of the 
isotropic Heisenberg ferromagnet, without energy conservation. A familiar problem one is 
faced with is the following. Whereas the upper critical dimension of this dynamic model 
is dj = 6, the one of its steady-state distribution, described by the usual |</>| 4 model with 
an n = 3 vector field <f>, is d* = 4. Thus the small parameter in which a dimensionality 
expansion can be made in the dynamic case is = 6 — d rather than = 4 — d, where 
d is the bulk dimension. For 4 < d < 6, the static critical behavior is given by mean- field 
theory and associated with the (then infrared-stable) Gaussian fixed point of the \4>\ 4 theory, 
even though the dynamic critical behavior is described by a nontrivial fixed point that is 
characterized by a non-zero value /* of the mode-coupling vertex and accessible to the 
expansion. 



3 



Unfortunately, this expansion is not tailored to capture the nontrivial static critical ex- 
ponents that emerge as d drops below 4. Therefore it is of somewhat limited use in the 
physically interesting three-dimensional case or, more generally, for d < 4. In order to 
to find out which scaling laws exist relating the critical exponents of dynamic bulk and 
surface quantities to known bulk and surface critical indices, it is essential to formulate 
the field-theoretic RG for fixed values of d < 4. We do this by extending existing massive 
RG schemes for semi-infinite systems^^ to dynamics. This yields RG (Callan-Symanzik) 
equations whose exploitation in conjunction with known boundary operator expansions!'!!® 
reveals that the dynamic bulk and surface critical exponents can be expressed completely 
in terms of known static ones, besides giving the scaling forms of quantities like the surface 
structure function. 

In order to check these finding we have performed high-precision computer simulations of 
a semi-infinite lattice model of classical Heisenberg spins whose dynamics is defined via the 
deterministic nondissipative equations of motion implied by their Poisson bracket relations. 
The advantage of this simple dynamics without noise is that recently developed extremely 
efficient spin dynamics alg orithmJlB can be employed to compute the temporal develop- 
ment from given initial spin configurations, which we choose from a thermal equilibrium 
distribution generated via a Monte Carlo simulation. 

It should be emphasized that this lattice model differs in an important aspect from the 
continuum model we consider: unlike the latter, it conserves the energy. Nevertheless, both 
models belong to the same universality class, as we intend to show. 

The remainder of this paper is organized as follows. In Sec. II both the semi-infinite 
lattice model (studied by simulations) as well as the semi-infinite extension of the continuum 
model J (utilized in our RG analysis) are introduced, and their dynamics specified. The 
definition of the continuum model also involves the specification of appropriate boundary 
conditions. We discuss this question first on a heuristic basis (Sec. |IIB 3|) . Going over to the 
path-integral formulation of this model in Sec. [HB4j , we then show in Sec. [II B 5| how the 
boundary conditions for both the order parameter (f> and the auxiliary (Martin-Siggia-Rose) 
field 4> can be justified in a systematic manner and derived from the boundary part of the 
dynamic action functional. Section |11 B 6| briefly recalls the fluctuation-dissipation theorem 
and discusses the meaning of some of the boundary conditions in this context. Section fTJ 
is devoted to the RG analysis of the continuum model. After giving the free response and 
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correlation propagators in Sec. |III A|, we explicate in Sec. [II IB] the renormalization of the 



theory, describe the massless renormalization scheme on which our subsequent RG analysis 
in 6 — e§ dimensions is based. To overcome the limitations of this scheme, we construct in 
Sec. |111 D| a massive RG scheme for fixed dimensions d with 2 < d < 4. The resulting Callan- 
Symanzik equations are given in Sec. E| and utilized to derive the scaling forms of the 



correlation and response functions. Details of our Monte-Carlo spin dynamics simulation are 
described in Sec. [TV]. Its results are presented and analyzed in Sec. [V[ Section |V| contains a 



brief summary and concluding remarks. Finally, in the Appendix arguments are given as to 
why the lattice model we study belongs to the universality class of our semi-infinite model 
J, even though it differs from the latter by conserving additionally energy. 

II. MODELS 

A. Semi-infinite lattice Heisenberg model 

The lattice model we consider is a classical isotropic Heisenberg ferromagnet on a d- 
dimensional simple cubic lattice whose sites i = (ii, ... , id), with i K = 0, . . . , L — 1 for 
K — 1, . . . , d, are occupied by spins Si = (Sf, a=l, 2, 3) of length \Si\ = 1. Free boundary 
conditions apply along the indirection and periodic ones along the remaining d—1 ones, so 
that the layers — and id = L — 1 are free surfaces. The Hamiltonian of the model reads 

Ti-ut = —J Si ■ Sj — Ji Si ■ Sj , (1) 

<i,j> (i,J> 

id or j d =£0,L-l i d =j d =o,L-l 

where the summations run over the specified sets of nearest-neighbor (nn) bonds 

The bulk and surface nn interaction constants J and J\ are ferromagnetic and measured in 

temperature units k&T. The dynamics is defined by the equations of motion 

dSi &H\a.t „ tri , 

which describe the precession of the spins Si in the local magnetic fields Hi oc —dH^t/dSi. 
They conserve both total spin ^\ Si (in the here assumed absence of external magnetic 
fields) as well as total energy E\ at = Hi a t[S(t)]. 

Conservation of magnetic energy is not normally considered a property of real ferromag- 
nets since the spin system can loose energy by processes not taken into account by Eqs. ([!]) 



and (0) such as coupling to phonons. In fact, in the continuum model J employed in our RG 
analysis, only the order parameter but not the energy is conserved. Arguments as to why 
both models represent nevertheless the same universality class are given in Appendix HJ. 

In our computer simulations, a d— 3 dimensional version of the above model is investi- 
gated. The equations of motion are numerically integrated for a given set of at least 
700 initial spin configurations generated by a Monte Carlo simulation of the thermal equi- 
librium state associated with Hamiltonian (^.000 Details of this simulation are explained 
in Sec. 

Quantities of primary importance for the interpretation of scattering experiments are the 
spin-spin cumulant 

C^(r-, Zl z';t-t') = (SUt)S^t>)) CUm 

= W)5?(0>- W)>(5?(0> (3) 

and its Fourier transform 

/POO 
d d-x re -i P r dte iujt C af3 (r;z,z';t-t') . (4) 
J — oo 

Here r = (ii— i[, . . . ,id-i— ^_i), while z = i 3 and z' = respectively. Further, t and i! are 
times to which the initial spin configuration at t — has evolved according to Eq. (Q). The 
average (.) is taken over the distribution of initial configurations. 

Specifically, we will be concerned with the dynamic surface structure function 

C^(p,u;) = C^(p;0,0;u;). (5) 

Before embarking on a discussion of its scaling properties and presenting our simulation 
results, it is useful to introduce first the continuum model on which our RG analysis is 
based. 

B. Semi-infinite model J 

1. Hamiltonian of the thermal equilibrium state 

The dynamic model we are going to consider is required to satisfy detailed balancelll 
and to ensure relaxation to a steady-state distribution corresponding to a thermal equilib- 



rium state oc e n ^ with the Hamiltonian 




Here the integrations extend over = {(x\\, z) G M. d \ z > 0}, the <i-dimensional half-space, 
and B, its (d — l)-dimensional boundary plane at z — 0, respectively. The order-parameter 
density <f> = (<f) a ) is a three-vector. 

Above d = 3 bulk dimensions, this static model is known to undergo at the bulk crit- 
ical point so-called ordinary, special, and extraordinary surface transitions, depending on 
whether the surface enhancement variable Cq is larger than, equal to, or less than a critical 
value c S p.0'i For d = 3, the surface cannot spontaneously order at the bulk critical temper- 
ature T c > because of the presumed continuous 0(3) symmetry of the Hamiltonian (^j). 
Hence only the ordinary transition remains in this case. Analogous statements apply to 
the lattice model ([!]), for whose d > 3 variant the role of the variable — Co is played by the 
'surface enhancement' (Ji/J) — (Ji/J) sp , where (Ji/J) sp is the critical value of the ratio 
J\j J pertaining to the special transition. 



2. Langevin equations 

Next we turn to the task of formulating an appropriate semi-infinite extension of the 
standard bulk model J. For reasons expounded elsewhere,B! we may assume that the surface- 
induced modifications of both the interactions as well as the dynamics are restricted to the 
immediate vicinity of the boundary B. Consequently, we use the stochastic bulk equation 

<j>(x,t)=\o(AH <p + foH <i> x < j>) + C(x,t) (7) 

for all points x with z > 0. Here £ is a Gaussian random force with average (£) = and 
variance 

(C(x, t) C'V, *')) = -2A 5 a ?A 6{x-x') 5{t-t') . (8) 

Further, Ti^ stands for the part of the functional derivative 

57~C 

— {x, t) = H^x, t) + 8{z) (c - d n ) d>(x, t) (9) 
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that remains away from the boundary plane B, namely 

H+(x,t)= (-A + r o + ^|0| 2 )</>. (10) 

The derivative d n in Eq. (|D is along the inner normal, i.e., d n = d z on B. 

In order to extend the model to the semi-infinite case, we must specify whether and how 
Eqs. (0) and (g) are to be modified in the vicinity of B. Owing to our locality assumption 
mentioned at the beginning of this sub-subsection, this should amount to a choice of bound- 
ary conditions for <p. For the sake of simplicity, we assume that the conservation of the 
order parameter is not violated by boundary contributions. This is physically reasonable 
since we took all bulk and surface terms of the Hamiltonian flf]) to have 0(3) symmetry, as 
is appropriate for a Heisenberg magnet whose interactions are isotropic even at the surface. 



3. Boundary conditions for <\> 

Building on previous work on model B§ii, we can now easily anticipate the proper 
boundary conditions. One boundary condition for <p should be the usual static one 

d n 4> = c <f>, (11) 

which ensures the vanishing of the contribution oc 5(z) of the functional derivative 

The other one is entailed by the required order-parameter conservation. This becomes 
clear if we rewrite Eq. (0) as a continuity equation: 

+ V- (j^+jf) = 0. (12) 

Here 

j(<*) = -A (VH r + f (^V W) (13) 

are the deterministic parts of the currents, and the noise parts satisfy 

C = -V-4 a) . (14) 

To ensure conservation of the total order parameter, no currents must leave the system. 
Hence the normal component of the currents should vanish, 

3 (») = n-j^=0, a = 1,2,3. (15) 



S 



If spin anisotropies were present at the surface (which is not uncommon), the conservation 
would be violated at the surface for some, if not all, components of 4>. 

Both boundary conditions, Eqs. (|TT| ) and (0), are valid in an operator sense, i.e., hold 
inside of averages over the initial values and the noise (yielding correlation and response 
functions). Note that the validity of the former has two immediate consequences: The surface 
contributions to the currents one would expect from the 5-function term of Eq. (Q) upon 
using the definition = — Ao V<57i/50 a rather than Eq. fll3|) disappears. Furthermore, 
substitution of the boundary condition (Jllf) into Eq. (|lf|) shows that the precession term's 
contributions (oc fo) to the currents j'"', a = 1,2,3, vanish, so that these latter boundary 
conditions become 

dji+ = d n (- A + t-o + y M 2 ) = 0. (16) 

The probability distribution of the noise clearly must also comply with the presumed 
order-parameter conservation. We prefer to discuss the consequences within the frame- 
work of the functional-integral (re-)formulation of the theory,§ilii where they manifest 
themselves as boundary conditions for the auxiliary or Martin-Siggia-RoseB (MSR) field 
introduced below. 



4- Functional- integral formulation 



The Langevin equations (0) can be rewritten as 



5(f> 

where 1Z = (7Z al3 ) denotes the reaction operator 



51~i 

cf>(x,t) = -['R,- — )(x,t) + C(x,t), (17) 



U aP = -A (5 af3 A + f <f) . (18) 



Since this operator acts on TC,p, which according to Eq. fllq) satisfies a Neumann boundary 
condition, the Laplacian it involves is self-adjoint on an appropriate space of (sufficiently 
smooth) functions satisfying this boundary condition. 

The measure e~ J ^V[(p, <f>] which appears in the equivalent functional-integral 
formulationcil of the theory can now easily be inferred. To this end, let us first recall 

which form the action J"[(j>, 4>] must have to ensure detailed balance and relaxation to the 



chosen equilibrium state. For the here considered case in which the noise has a Gaussian 
probability distribution, this isEEEl 



J 




(19) 



where a pre-point discretization of time is understood to be employed. 

The action of the bulk model J is known to be of this form, with the reaction operator 
1Z being given by Eq. (|T8|). If we accept the boundary conditions ( |TTD and (|16D, then 
contributions to the action that are localized on the surface vanish. Consequently this result 
for the action must also hold in the semi-infinite case we considered, with J x interpreted as 
the volume integral f Rd . 

Conversely, one can start from an action of form ( |I9"D and derive the boundary conditions 
in a systematic fashionil along the lines followed in Refs. ^ and ^. The various general 
assumptions we have made (consistency with bulk model J, only local modifications of the 
dynamics at the surface, absence of nonconservative surface terms, etc.) can be combined 
with relevance/irrelevance considerations to conclude that the reaction operator reads 

K = A ( S a/3 V V - fo e a ^ (j>A , (20) 



where V acts to the left. This is identical to Eq. fll8 ) up to the replacement of the Laplacian 
by the symmetric expression VV. 

The substitution of this form of 1Z into Eq. (|l^) and an integration by parts (making no 
use of the boundary conditions ( |TTD and (|16|)) yields 

J = [°°dt( f (0-(</)-Ao/oW x0) 

J-oo \JRf 1 

-A (H^-tyA^} 
-A o ^{(A0) .(co-d n )<l>+ 

+ 5(z) (c - d n )<f> - f x $\ ■ d n ^j . (21) 

The singular piece oc S(z=0) present in the boundary integral J B is familiar from Refs. |7]|||9|- 
It results from the coincidence of two 5 functions. This singularity can be cured by replacing 
one of the 5 functions by a smeared-out smooth analog such as Sg(z) = B e~ Bz , with a large 
positive but finite value of B. 
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5. Boundary conditions as boundary equations of motion 



Starting from the action fll9|), we can now obtain the boundary conditions for both (f) 
and in a standard mannerB! as 'boundary contributions to the equations of motion'. This 
works as follows. We add source terms Jj [<fi, </>] to the action and consider the generating 
functional Z[j] = J'T>[<f), <fi] e~^ - ^, where j stands for the set of all sources considered, 
including eventual ones localized on the surface. That is, the source part of the action can 
be written as 



JA4>, <t>] 



dt 



j K o K + / i P o 



(22) 



where K = K (x,t) and O p = O p (x\\,t) are local functionals of and </>. From the 
invariance of the generating functional Z[j\ with respect to a change of variables <p — > 4>+5<fi 
and (p — > (ft + 5(f) with arbitrary (smooth) functions 5(f) and 5(f) we may then conclude that 



(5 J + SJ J ) J = 



(23) 



Here 5 J and bjj are the implied changes of first order in 5<p and 8(f) of the action J and 
its source part J}, respectively, and denotes the average in the presence of sources. 
Explicitly, we have 



8J 




J(p5(j) + J^54> 

J<p B 5<t> + Jj, B 5<t> 
+ Jd n <p d n 5(j) + J 9n ^ dJ4> 

+ J(A4>) B + ^(A0) s A6( t> 



(24) 



with 



U 4> 



-d t (f) - A 



A + C/o 



A(f) 



— fo 



AUx (j) +0 x H 



Jj, = d t ^-X A(n 4f -2^ -Ao/o^x (j) 

J*b = (c -d n )^-\ [U 2 + c 5(0)}d n 4) 
+X fo4> x (co - #„)0 , 



(25) 



(26) 



(27) 
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Jj B = A d n (2(f) - HA + A /o 4> x d n cf) , 



(28) 



Ja n 4, = *o6(0)d n $, (29) 

J dn 4> = ^m(d n -c )4>, (30) 

*7(a</») b = A o <9 n (31) 

and 

^ (A0) B = A o {d n - c o )0 • (32) 

The subscript B at local quantities indicates their restriction to the boundary. Further, U 2 
is the matrix of second derivatives of the \4>\ 4 potential: 

Uf = r 5 al3 + ^ (5 al3 + 2 V) , (33) 
6 

and for convenience, we have introduced the field 

$ = 7?, T -0 = Ao(-A0 + /o0x0) , (34) 

where 1Z 1 denotes the transposed reaction operator, i.e., (1Z T ) al3 = 1ZP a . We shall refer to 
$ as response field since this is its known physical sig nificanceBi (as can be read off again 
from our Eqs. (|41~D and (|42|) below); it should not be confused with the auxiliary field <fi. 

For simplicity, let us include here in the source part of the action merely bulk sources 
J(x,t) and J(x,t) as well as surface sources Ji(x\\,t) and Ji(x\\,t) which couple to <f>, <f>, 
0B) and <pBi respectively. Owing to the arbitrariness of 5(f) and 5cf>, it follows from the above 
results that the 'equations of motion' 

J<t,{x,t) = J(x,t) , x£B, (35) 

Jj{x,t) = J{x,t) , x£B, (36) 
and the 'boundary equations of motion' 

J p (x\\f,t) = J 1 (x ll ,t)S P;<Pl3 + J x {x h t)5 pt fa , 

(37) 

P = <t>B, 4>b, d n </>, d n 4>, (A<p) B , (A0) B , 
12 



hold inside of averages with the action J + Jj. From the latter and Eqs. (|27Q-(|32D , we get 
the previously given boundary conditions for <f>, Eqs. ([EL]) and (|T6|), and two additional ones 
for <fi: namely 



and 



d n 4> = (38) 
(d n -c Q )$ = 

{d n - c ) A (- A0 + f 4> x </>) = . (39) 



6. Fluctuation- dissipation theorem 

The significance of the boundary conditions ([3£]) and (|39|) can be understood as follows. 
Let us add a (possibly time- dependent) magnetic field term to the Hamiltonian making 
the replacement 



n^n h = n- h-4>, (40) 

where we assume h(x,t) to have support only off the surface. This induces the change 

J^J h = J-! h $ (41) 



of the dynamic action. Hence we recover the usual correspondence (known from the bulk 
casell) 

*(br < 42 » 

between functional derivatives with respect to h(x), taken at h = 0, and insertions of the 
response operator on the right-hand side. 

Furthermore, the fluctuation-dissipation relation 

-9(t - t') (r(x, t) t')) = (4>"(x, t) $ v, 0> (43) 

can be derived as in Refs. [H], |22|, 24, and 32, owing to the form (19) of the action. 



The significance of the boundary condition (|39|) becomes clear if we let x' in Eq. (|43| 
approach a point on the boundary plane B: it ensures the consistency of the fluctuation- 



dissipation theorem (|I3|) with the boundary condition (11 
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To understand the Neumann boundary condition (j38| ), note first that according to 
Eq. (|i9|), the reaction operator couples V0 to the current operator — X^VSTi/ 5cf). This 
boundary condition ensures that the reaction operator is self-adjoint. In addition, it can be 
combined with Eq. ( [ZED to see that the boundary equation of motion Q57D for p = (fig leads 
back to the boundary conditions ( |TED for the currents. 



III. RG ANALYSIS OF THE SEMI-FINITE MODEL J 



A. Preliminaries 



We now turn to the RG analysis of the semi-infinite model J introduced in the previous 
section. To this end, two RG schemes will be used: a massless one based on minimal 
subtraction of poles and the expansion about six dimensions, called RSi, and a massive one 
for fixed dimensions 2 < d < 4, called RS2. 

Before embarking on a discussion of either one of these, we must set up some notation. 
Let us define the generating functionals of (connected) correlation and response functions 

W[J, J,K,Jt, JJ 

= In ^ e (j>&+( J >&+( K ^ x <P)+( J i>4>B)+(Ju<pB)^ (44) 

and 



g[j, J;Ji, Ji] 

= In < / e (^*)+(- 7 ^)+( J i.*8)+(Ji,08)^ ? (45) 
where we have introduced the convenient short-hands 

/OO /' 
dt / d d xJ{x,t)-(f)(x,t) (46) 
-00 Jw 1 , 



and 

"OO 



/OO /' 
dt / d d - x x\\ J x {x h t) -<t>B{x h t) . (47) 
-00 JB 

For the cumulants generated by these functionals we write 



W (N,N,L;M,M) 



N N L M M run: 

n n ^ x ^ i n * n * 

j=l k=l 1=1 m=l n=l 
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(48) 



and 

q(N,N;M,M) 

I N N M M v cum 

n^n^n^n*) (49) 

,7=1 k=l m=l n=l ' 

We normally suppress the tensorial indices &i, . . . , /3m of these functions, but give their 
space and time (or momentum and frequency) coordinates (suppressed above) when dealing 
with specific ones. 

The reason for considering the functions w( N ' N ' L ' M ' M ^ should be clear: aside from multi- 
point cumulants of the basic fields <f> an d </>, insertions of the composite operator cp x cf> are 
needed because it appears in the fluctuation-dissipation relation (|43|). 

The free response propagator and free correlator one needs to compute the W and G 
functions defined above are the same as for model Bb, and for the case To > 0, may be 
gleaned from any of the Refs. [F|,f§,^| In a mixed pzu representation (where p G M.^ 1 is 
a (d — l)-dimensional wave- vector conjugate to xn while u denotes the Fourier frequency 
variable associated with t), the free response propagator reads 

- - — [ e ~ K +\ z -' z \ - f + e -K+iz+z) _ g+ e -(*+*+*-i)] I ( 5 Q) 

with 

f _ f , K . k) = K ± k t (4 - 4) ~ c ° t K ± ( k2 ~ 4) + « T - 4)1 (51) 



and 



= ± ( K± « • Co «) = 2c k± (/t 2 - 4) 

T ' ' K± (k 2 - 4)( c + k t) ~~ K T ~ 4)( C + K ±) 

Here 

« = + vV + r o , (53) 
and k± denote the roots with positive real parts of the equation 

4 = P 2 + (ro/2) ± [(ro/2) 2 + i (cj/Xo)] 1/2 . (54) 
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The free correlator C 5 a/3 = (0 a 0^) cum can be expressed in terms of the free response prop- 
agator as 

POO 

C \p; z , z' '; lo) = 2A / dz G(p; z, z;u) 
Jo 

x(p 2 -d!)G(-p; z',~z; -u) . 

(55) 

Owing to the presumed 0(3) symmetry of the Hamiltonians ([!]) and (|J) of our lattice 
and continuum models, the only surface transition that is possible in three dimensions is 
the ordinary transition. We can benefit from the fact that its asymptotic critical behavior 
can be studied by taking the limit Cq — > oo (see, e.g., Ref. |5] and below). The simplified 
expressions for the free response propagator and correlator which then apply correspond to 
the replacements of the coefficients f± and g± by 

/± = lim /±(k±, k t ; c , k) = _1±_^± ( 56 ) 



and 



2k± 

= lim (/±(k±, k t ; c , k) = ■ , (57) 



respectively. 



B. Massless renormalization scheme near six dimensions (RSi) 

We here restrict ourselves to bulk dimensions 4 < d < 6. Then the static critical behavior 
is described by Landau theory. The Gaussian fixed point, uq = 0, of the |0| 4 theory is 
infrared-stable. In part of the calculation uq therefore can be set to zero. This is possible 
as long as we consider quantities that have a nonsingular and nonvanishing limit u$ — > 0. 
However, we must keep in mind that the linear scaling field u associated with uq is dangerous 
irrelevant .IS Quantities like the free energy density or the spontaneous magnetization vary 
as inverse powers of uq ~ u for uq —>■ 0, and hyperscaling is broken. Accordingly, already a 
full scaling description of the static bulk critical behavior requires the inclusion of a second, 
so-called thermodynamic, length besides the bulk correlation length. Finally, in applying 
techniques of renormalized field theory, we must remember that both the static as well as 
the dynamic theories are not renormalizable for d > 4 if uq ^ 0. Single insertions of the 
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local operator to which Uq couples can be renormalized, but the additional counterterms to 
which it gives rise are not sufficient for curing the additional ultraviolet (uv) singularities 
produced by multiple insertions. 

We now consider the W functions (EST) , where we restrict the temperature T to values 



above the critical temperature T c and set Uq = temporarily. Since the Hamiltonian (^) 
becomes Gaussian for uq = 0, there are no static uv singularities to cure. Hence no amplitude 
renormalization of the order parameter is required, i.e., Z^ — l. By power counting one finds 
that counterterms of the form J dt j Rd (f> ■ (p are also not needed. Since this implies that 
the product of amplitude renormalization factors for (f> and (f> is unity, an amplitude 

renormalization of is not required either {Zi = 1). A (bulk) counterterm oc <fi ■ A(j) is 
ruled out for u = 0. If u did not vanish, the O(u ) contribution to such a counterterm 
would diverge as A 4 for d = 6 within a theory regularized by a cutoff A. 

We assert that the following counterterms are sufficient to renormalize the generating 
functional W: aside from those implied by the reparametrizations 

\ = l x- 4 Z x (f,d)\ (58) 

and 

fo = H S ?[ZxV*,dj\- l f, (59) 

only a counterterm of the form J dt J Rd K-A<f> is required, where /i is an arbitrary momentum 
scale. More precisely, we claim that the cumulants generated by the functional 



/oo r 
dt K-A4> (60) 
-oo JR^ 

are uv finite when expressed in terms of A, / (and Cq or its dimensionless equivalent c = Cq///). 

This conclusion is based on the following observations. The detailed-balance form of the 
action (|19|) in conjunction with the constraints imposed by the conservation of the order 
parameter and power counting restricts the possible bulk counterterms to those included 



in Eq. (60). Using this result as input, one can consider the renormalization of the bulk 



analog of the fluctuation-dissipation relation (fig) . For convenience we employ dimensional 



regularization and fix the counterterms by minimal subtraction of poles at d — 6. From 
the uv finiteness of the correlation function of the renormalized function on the left-hand 
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side of this equation two conclusions may then be drawn: the renormalization factors of fo 
and Ao are reciprocal to each other, and the renormalization function of the -RT-dependent 
counterterm is related to Z\ in the stated fashion.0 The result means that insertions of the 
response field renormalize just as (p, requiring no additional counterterms. 

A final step remains to complete the argument: we must show that no further surface 
counterterms are needed. Given the causal structure of the theory (according to which 
at least one (f) must occur in any monomial of the action), power counting restricts the 
possible candidates for such counterterms in d = 6 dimensions to boundary contributions 
to the action involving monomials of the form Ao </> 2 , Ao <fi ■ <j>, and Ao0 a <^</> 7 (as well as 
similar ones with derivatives), where the coefficients have momentum dimensions 1, 3, and 
0, respectively. Now the cubic nonlinearity of the bare action can be rewritten as 



oo 



oo 



J mcv = Ao fo e a ^ / dt / (Vr) / (61) 



upon integrating by parts and utilizing the boundary condition (|TTD. Thus each leg of the 
vertex oc fo comes with a derivative V. This reduces the superficial degree of divergence of 
such uv boundary singularities with two <fi legs by two, making it negative (uv superficially 
convergent). By a similar argument, surface counterterms involving one <fi and two <fi fields 
are ruled out. Hence we are left with surface counterterms oc • and analogous ones 
with up to two additional derivatives. 

To proceed we follow Ref. |22| and perform the integral dt of the fluctuation-dissipation 
relation ( |43f) . This yields 



{Wx)^)* = -A o A'<0 Q (a ; )^(a ; ')>_ o + Ao/o<0 Q (x)(0x0)^xO> Lj=o , (62) 

where the superscript 'st' indicates a static quantity while the subscript u at the expectation 
values on the right-hand side means their Fourier transform with respect to the time differ- 
ence t — t' . We multiply this equation from the right by the inverse of the static propagator 
on the right-hand side and from the left with the vertex function Tia^. The result is@ 

(x, x'; uj=0) = Ao (- A' + r Q ) [ - A 6(x - x') + f 1> ; { ^ x0) , (x, x'; uj=0)] , (63) 

where F^.^^p means a vertex functions with a single 4> a leg and an insertion of the com- 
posite operator (<fi x (p) 13 . Owing to the operator — A' + r (inverse static propagator) on the 
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right-hand side, the uv behavior of possible primitive local divergences of the vertex func- 
tion on the left-hand side is improved by two powers of A. This is sufficient to ensure that 
no additional surface counterterms with one and one leg are needed. Since the vertex 
function r^ Q .^ x ^ /3 has an explicit V acting on the external leg, it has a primitive loga- 
rithmic bulk singularity that is cured by the subtraction provided by the bulk counterterm 
oc K ■ A<p. On dimensional grounds one might anticipate logarithmically divergent surface 
counterterms of the form (d n <p) ■ (tq — A)0 and </> • <9„(ro — A)<p, but these are annihilated 
by the boundary conditions (|38l) and (|16|), respectively. Note also that the restriction to 
uj = is unproblematic here because each factor of u (i.e., each time derivative) reduces the 
superficial degree of divergence by four. 

A further comment is appropriate here. Field theories for systems with boundaries are 
known to have the following feature: Besides one-particle irreducible (1PI) graphs, one- 
particle reducible (1PR) ones may also require 'final subtractions' and hence contribute to 
renormalization functions (cf. Sec. II. B. 6 of Ref. |5|). Nevertheless our above reasoning based 
on 1PI graphs is conclusive since the power counting would not be changed if we contracted 
1PR graphs whose external free legs are amputated to a point. Hence the counterterms 



included in W rcn are indeed sufficient to cure the uv singularities of the W functions (pEBD- 
The uv finiteness of the G functions generated by 

G ien [J, J, K, J 1; Jj = g [J, J, K, J 1; Jj (64) 

when expressed in terms of A and / (as well as tq and cq or their dimensionless equivalents 
t = r /fi 2 and c = Co/ fx) follows as a simple corollary from the fact that the response field 
<& renormalizes just as <p (namely trivially). 

C. RG analysis in 6 — e dimensions 

Utilizing the results of the previous subsection, one can perform a RG analysis of quan- 
tities that are finite and nonzero for uq = 0. This criterion is satisfied by both the W and 
the G functions for r > 0, as can be checked via perturbation theory in f . Since such a 
RG analysis is standard we can be brief and just state its principal results. 

To one-loop order the renormalization function Z\ is given byii 

z ^- l -^& i + 0(fl) - (65) 
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Equation fl59|) implies that the beta function 

W) = ^\ f (66) 
can be written in terms of the exponent function 

77 A (/) = /^| lnZ A , (67) 

namely 



W) = -f 



6-d , ' 
■ w) 



(68) 



2 

where d fl \ denotes the derivative at fixed bare coupling constants fo, tq, and cq. The 
infrared- stable fixed point /* for 4 < d < 6 is given by the nontrivial root of the equation 
/?/(/*) = 0. From Eq. (|6~8"D we find for the value of the exponent function r]\ at the fixed 
point /* the result 

vl = Vx(f*) = ^, (69) 
which we insert into the general expression 

3 = 4-7^ (70) 

for the dynamic critical exponent 3. Since the correlation exponent 77 is zero for d > 4, the 
final result 3 = (d + 2)/2 for 4 < d < 6 is consistent with the established valueSSIIlIIi] 

3 = ^1^, 2<rf<6. (71) 

The latter result is known to follow most easily from the observation that the characteristic 
frequency of isotropic ferromagnets for T > T c is the Larmor frequency and hence scales as 
the (static) bulk magnetic field 

Next we turn to the analysis of the critical behavior of surface quantities. We restrict the 
following discussion to the case of the ordinary transition, the only surface transition at bulk 
criticality that remains in the physically interesting three-dimensional case. The surface 
enhancement variable c= c^/fi transforms as c —>• c(l) — c/i under scale transformations 
fi — > fi£ and hence approaches the fixed-point value c* rd = 00 in the infrared limit I — > 0, 
provided its initial value is positive. In this case we can set Cq = 00 from the outset. Surface 
quantities involving the surface fields </>g or #g then vanish. 
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Let F ren be any of the two renormalized functions Wi^' N ' L ' M ' M ^ and Grm' 7V,M ' M ' 1 generated 
by the functional (|60|) and (|64l) for uq = and r = tq/^j 2 > 0. The invariance of the 
regularized bare functions with respect to a change of /i implies the RG equation 

+ p f d f - 2rd T + (4 - r]\)\d\ - c«9 c ] F ren = . (72) 

This may be utilized in a familiar manner to obtain the asymptotic scaling forms of these 
functions. The result one obtains for the pair correlation function C a/3 = S a)3 C (= G^ 0,2 ' 0,0 )) 
at T = T c agrees with the more general one predicted in our previous paperjl! 

C(r; z, z'- 1) ss r 2 - d ^ T (z/r, z'/r; tr~*) , (73) 

if the classical value r\ = and the implied one 3 = (<i+2)/2 are substituted (as is appropriate 
for A < d < 6). Here we have suppressed the variables \i and A, setting both to unity. The 
variable c does not appear on the right-hand side because the scaling function T is a property 
of the c = 00 fixed point. Deviations of c from the value c = 00 produce corrections to the 
displayed leading infrared contribution. 

The derivation of the scaling form of the surface correlation function 

C u (r;t)=C(r;0,0;t) (74) 

is not quite so straightforward because C\ c=oc vanishes whenever z = or z' = 0, as a 
consequence of the Dirichlet boundary condition into which Eq. ( |Tl"|) turns for cq = 00. 
One possibility to cope with this difficulty is to the study behavior of C\\ in the limit 
c — > 00. As is expounded elsewhere,B0&i@ this can be achieved by an expansion in powers 
of l/c . According to the boundary conditions (|Tl"D and (p9|), the boundary operators cf>B 
and <E»b can be replaced by Cq l d n (f>s and Cq respectively. From previous detailed 

investigationsi'iS of the 1/cq expansion it is therefore clear how the scaling forms that the 
correlation and response functions involving these boundary operators take at the ordinary 
transition can be determined: making the replacements (ps — * d n (f>B and $g — > d n <&B, one 
studies the so-obtained analogs of these functions, with cq set to 00. 

An alternative strategy, which leads to equivalent results, is to use the boundary operator 
expansion (BOE). According to Eqs. flTTD and (0), both the order parameter density <fi as 
well as the response field $ satisfy Dirichlet boundary conditions if Cq = 00. In analogy 
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with the static case, a BOE of the form 

4>(x, t) « D(z, t) d n (f>(x B , t) + ... (75) 

Z— >0 

and its counterpart involving $ and d n $> should hold clS X — ^£C|| ; Z ) approaches the surface 
point xs = (aj||,0). (We have suppressed contributions proportional to the operator 1, 
which occur when the average (0) does not vanishBtiH) Equation [7^ implies that cumulants 
involving the surface operators d n (p and <9 n $ give access to surface correlation functions. 
We refrain from doing this within the framework of the €q expansion and turn directly to 
RS 2 , the massive RG scheme. 



D. Massive renormalization scheme for 2 < d < 4 (RS2) 

Our aim here is to extend the massive RG scheme for semi-infinite systems developed by 
one of us and ShpotH'il to the dynamic theory of the model- J action fl2~l"|). We assume that 
2 < d < 4 and give up the restriction w = 0, i.e., both u and f are assumed to be nonzero. 

The scheme can be extended for general values of cq. However, as we are primarily 
interested in the dynamic surface critical behavior at the ordinary transition in d = 3 bulk 
dimensions, we can simplify the analysis by setting Co = 00 from the outset. The advantage 
of doing this is considerable: for general values of Cq, the renormalization factors associated 
with surface operators (called 'surface renormalization factors' for short) depend on the 
renormalized coupling constant u (to be defined below) and the ratio c/m, where c and 
m are the renormalized analog^! of the bare surface enhancement Co and the renormalized 
mass m (to be introduced below), respectively. This additional dependence on the mass ratio 
c/m makes the RG analysis rather cumbersome. If we set Co = 00, we focus directly on the 
asymptotic regime c/m = 00 and avoid this difficulty because the surface renormalization 
factors (which are of purely static origin) can then be chosen to depend merely on u. 



1. Static bulk renormalization functions 

Let r^ 7 -* be the static bulk vertex function with iV legs of type <ft and I insertions of </> 2 /2, 
and f (<7, Q) the Fourier transform of this (translationally invariant) function, up to the 
momentum-conserving factor (2ir) d 5(^2 q + ^2Q)- Here q and Q are the N and I momenta 
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conjugate to the positions x and X of the legs and the inserted operators, respectively. We 
write 



t = m 2 + 5m 2 , 



(76) 



<(> = [Z^u)] 1 / 2 c/> VCQ , 



(77) 



2 = [Z^(«)]- 1 [0 2 ] rcn , 



(78) 



and 



u = Z u {u) m u , 



(79) 



which introduce the renormalized mass m, the renormalized densities (f) Ten and [</> 2 ] ren , and 
the renormalized coupling constant u. 

The mass shift 5m 2 and the renormalization ( l Z J ) factors are fixed via the familiar nor- 
malization conditions 



q=0 



m 2 , 



(80) 



^ H 2 ) i \ 



1, 



9 =0 



(81) 



q=Q=0 



and 



f£b,ren({*};W^' 



{<Z;=0} 



for the renormalized static bulk (b) vertex functions 

rSU-,rn,u) = [Z^u)] N/2 [Z^ (u)] 1 T^\-; r , u ) 



(82) 



(83) 



(84) 



with (N, I)± (0,1), (0,2). 
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2. Dynamic bulk renormalization functions 



We introduce renormalized auxiliary and response fields (f> ren and 3> ren , a dynamic bulk 
renormalization factor Z\, and the renormalized bulk variables A and / such that 

<(> = [Z^u)}- 1 / 2 4> Ien , (85) 

*=[^(ii)] 1 /2* renj (86) 

X = m~ 4 Z x (u,f)X, (87) 

and 

f = m^[Z x (uJ)]- 1 [Z^u)] 1 / 2 f. (88) 

To fix the function Z\, we choose the normalization condition 

d 



= -2m" 4 A5 a/3 , (89) 

g=oj=0 



where f ^'^ e /3 n ' ) denotes a renormalized dynamic bulk vertex function in the momentum- 
frequency (qui) representation. 

Let us add a few clarifying remarks. Note, first, that the renormalization functions 
introduced above are sufficient to absorb the uv singularities of the vertex functions with 
arbitrary numbers N and iV of and legs of the dynamic bulk theory for d < 4. Hence 
the above reparametrizations also yield uv finite renormalized functions when applied to 
the bulk analogs of the N + N point cumulants with iV cf) and iV <p fields , i.e., to the bulk 
analogs of the functions W^ N ' N ' ' ' ' defined in Eq. (|48|). The same remark applies to the 
bulk analogs of the G functions (fE3). 

The meaning of the multiplicative renormalizability of the response field $ with respect 
to the renormalization of the N + N point bulk vertex functions with an insertion of the 
composite operator 0x0 has been discussed elsewhere!l@ and need not be reiterated 
here: it implies, in particular, that the renormalization of ~ ^ involves a subtraction 
proportional to chosen in conformity with the multiplicative renormalization d86|) . 



Our second remark concerns the renormalization of ren and / in Eqs. (85) and (j38|). 



The fact that no primitive uv singularities involving a counterterm proportional to 

24 



~ — 1/2 

occur implies that the renormalization factor of <p rcn is given by Z, (up to a uv finite 
factor). This result means that the latter product of operators transforms according to its 
engineering dimension under RG transformations and is related to the conservation of the 
order parameter .0 

The form of the renormalization factor of /, Z^ 1 z\ , follows from the fluctuation- 
dissipation theorem (]4"5|); it tells us that (p a /Q a is a RG invariant. As a direct consequence, 
a relation generalizing Eq. (^) holds between the beta function f3f and exponent functions, 
which now are defined via 

p f (u, f) = md m \ f , p u (u) = md m \ Q u (90) 

and 

Vk{u, f) = md m \ Q \nZ K , k = A, 0, 2 , (91) 



respectively. We have 



P f (u,f) = -f 



6 - d , . rj4>{u) 
— 2 V\(u,f) + ^- 



(92) 



Since rj = r)<f,(u*), where u* is the nontrivial zero of f3 u , this form ensures that the established 
result (|7I| ) for the dynamic exponent 3 is obtained if the value 

^=^r) = ^ (93) 

pertaining to the infrared-stable fixed point (u*, /*) is substituted into Eq. fl70f) . 

Finally, let us note that the renormalization factors Z^, ^2, Z u , Z\ introduced above are 
all uv finite for d < 4, although they are logarithmically divergent in the uv if d = 4 (i.e., 
they have pole terms at d = 4). In other words, if ci < 4, then the uv singularities of the 
(static and dynamic) bulk theory are absorbed by the mass shift 8m 2 . 



3. Static surface renormalization functions 

In order to generalize the above approach (RS2) to the semi-infinite case, we set Cq = 00 
and consider the analogs of the G functions Q4"9| ) that result from these when the boundary 
operators and </>g are replaced by the normal derivatives d n <& and d n (f>, respectively. We 
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denote these functions as 

q(N,N;M,M) 



I N N M M 

= ( n ^ n ^ n n ^ 



k=l 



m=l 



n=l 



(94) 



CQ =00 



Following Refs. |25| and |26|, we introduce the static surface renormalization factor Zi >ao (u) 
and the renormalized surface operator (d n (f)) Tcn via 

d n <f> = [Z^u) Z ljOQ (u)] 1/2 {d n <t>) ren . (95) 

Next we take the normal derivative on both sides of the fluctuation-dissipation relation 
( fi3|) with respect to x' and set t' — 0. The result reads 

-9(t) (<P a (x, t) d n (^(x B , 0)) = {r(x, t) d n ¥(x B , 0)) , (96) 

which suggests to renormalize <9 n 3> in complete analogy with Eq. ( p5|) by 

d n $> = [Z^u) Z hoo (u)} 1/2 (d n *) ven . (97) 

This definition ensures that the modified fluctuation-dissipation relation (|9~6"D carries over 
to the renormalized theory. Moreover, it establishes consistency with the renormalization 
of the corresponding static correlation function in Ref. provided we fix Z l oo (u) as in 
Eqs. (7.10a-10b) of that reference. To this end we define the renormalized functions via 

N+N+M+M M+M 

GgMM) = Z~—^ Z~— GgMM) , (98) 
if (N, N; M, M) + (0, 0; 1, 1). The excluded function 



(99) 



requires an additive counterterm, which we choose such that the normalization condition 



p=LU = 



holds. To specify Zi yOD , we require that 

d 



_bLG(0Ai,i)fp u ) 



p=uj=0 



1 

2m 



(100) 



101) 



Equations ( |100|) and (101) are equivalent to the normalization conditions (7.10a,b) of 
Ref. |2§ together with Eq. (g9[) they imply that 

requires a subtraction and the 
renormalization factor Z\ ^ is the same as in the static case. 
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E. Callan-Symanzik equations 



The Callan-Symanzik (CS) equations can now be derived in a standard fashion. We take 
a derivative d TQ of the bare (dimensionally regularized) functions (P^) at fixed values 
of the other bare interaction constants Uq and /q. Using the above reparametrizations and 
definitions of the beta and exponent functions then yields 



m N + N + M + M M + M 

V m + 7fo + T]i : 



/~i(N,N;M,M) _ r>(N,N;M,M) (Tr\9\ 



with 



= ™^ + ^Yu + p t§f ~ (4 ~ Vx) xdx (103) 

and 

<r™ M) = ^^^^-J +M)/2 (m 9 m | r ) a T0 G^ W; ^ M) 

= (2-^)m 2 [9 T0 GL™ M) ] ren . (104) 

Here the exponent function r)i t00 (u) is defined by setting k = (l,oo) in Eq. (PT|). Just like 
the other static functions r) u (u), rj^u), 77^2 (it), it depends only on u (and d), but not on the 
dynamic coupling constant /, and is precisely the same as in Ref. The inhomogeneities 
Roo,'ren M ' MS> involve renormalized functions with an insertion of — J d d x<j) 2 /2, given by 



We proceed along lines similar to those followed, for example, in Refs. [50] and [51], in order 
to derive the asymptotic scaling forms of the response and correlation functions from the 
CS equations (|102| ). The deviation 5tq of the bare variable tq from its bulk critical value tq c 
depends on the temperature difference r = (T — T c )/T c according to 

St = T - T 0c ~ T , (106) 

which holds if r is sufficiently small. Near criticality the mass m — i.e., the inverse £ _1 of 
the (second-moment) bulk correlation length £ — behaves as 

m ~ I r - r 0c r with v = (2 + r^) -1 • (107) 
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Furthermore, integration of Eqs. (|91|) and ( |9~0"D gives the asymptotic dependence 

Z K ~ m ^, k = (ft, (p 2 , A, (1, oo) , (108) 



for (u, f) — > (u*,f*) or m — > 0. We insert this result for k = A into Eq. (p7[) , substitute 
expression d93|) for 77^, and arrive at 

A~m 3 A , (109) 



where 3 is the dynamic bulk exponent (|7Tj). With the aid of these results, it is straightforward 
to deduce the scaling forms 

g£> n ^ m \{x}^t 07 u Jo,\o) 

(110) 

Here (3 = {y/2)(d— 2+77) is a standard bulk critical index, while /3° rd is its surface counterpart 
for the ordinary transition. (For recent estimates of the numerical value of the latter at d = 3, 
see Ref. its references, and Ref. [31].) The set {x} comprises all position coordinates on 



which the respective function depends. The case (N, N; M, M) = (0, 0; 1, 1) is special in 
that the term 8(t— t) S(xn— x») G^ 0,0;1 ' 1 - l (p=0, uj=0), which results from the subtraction in 
Eq. (|9~9"D, should be subtracted on the left-hand side of Eq. (|110|) . We have suppressed this 



term, because we consider G^ ' ' 1 ' 1 ) here not as a distribution, but as a function for t — t > 0. 

Let us choose (N, N; M, M) = (0, 2; 0, 0) in Eq. ( |110| ) and consider the case of the spin- 
spin cumulant (Q). If we set Ao = 1 for convenience, we obtain the scaling form given by 
Eq. ( |T3"D in the limit m — > 0. 

The scaling form of the surface structure function ( [7^D at the ordinary transition can 
be derived from the expansion 

of G (0,0;0,2) 

in powers of c 1 (see Eq. ( |110|) ). Alternatively, 
we can combine the CS equations ( |102|) with Eq. ( |110| ) and the BOE ( |T5f) (applied to the 
renormalized theory) to conclude that the coefficient function D(z, t) asymptotically satisfies 
the CS equation 

771,00 



T>r. 



D(z,t) = 0. (Ill) 



2 

In the limit m — > 0, Eq. ( |1 1 1| ) yields a leading short-distance singularity of the form 

D(z,t) w D z 1+ <~l 2 , (112) 
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This in turn implies that the Fourier transformed surface structure function Cn(p,u) at T c 



can be written as 



C n (p,uj) «p'' 




(up } ) . 



(114) 



The limit p — > exists. By consistency, we must therefore have 



Cii(0, uj) = const uj 



( 3 +l— r?f d )/3 



(115) 



In the next section, we will check the predictions ( |114j ) and ( |115| ) by means of accurate 



Monte-Carlo spin dynamics simulations. 

IV. MONTE-CARLO SPIN DYNAMICS SIMULATION 

The Monte-Carlo spin dynamics simulation works as follows: A Monte-Carlo simulation 
of the lattice model with the Hamiltonian ([!]) yields a spin configuration that is used as 
initial condition for the integration of the equations of motion (0). When the integration 
is completed, the time evolution of the spin configuration is analyzed for position and time 
displaced correlations. The correlation functions are then stored in arrays, and a new initial 
condition is generated by the Monte-Carlo simulation. Typically, this is repeated 700 to 
1000 times. The correlation function C a ^(r; z, z'; \t—t'\) [cf. Eq. (||)] is finally obtained by 
averaging over the individual measurements. 

The Monte-Carlo algorithm is chosen as a hybrid scheme consisting of Metropolis sweeps, 
Wolff single cluster updates,il and overrelaxation sweeps.il Typically, one hybrid Monte- 
Carlo step involves 10 individual steps, each of which can be one of the updates listed 
above. Both the Metropolis and the Wolff algorithm work in the standard way, where the 
reduced coordination number of the lattice at the surfaces and the modified surface coupling 
J\ must be taken into account. The acceptance probability p of a proposed spin flip in the 
Metropolis algorithm is defined by p(AE) = l/[exp(AE/kBT) + l], where ks is Boltzmann's 
constant and AE is the change in configurational energy of the proposed move. 
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The overrelaxation part of the algorithm performs a microcanonical update of the con- 
figuration by sequentially rotating each spin in the lattice such that its contribution to the 
energy of the whole configuration remains constant. The implementation of this update 
scheme is straightforward. To see this, note that the energy of a spin with respect to its 
neighborhood has the functional form of a scalar product according to Eq. ([!]). There- 
fore, a spin can be rotated about the direction of the local field generated by its neighbors 
without changing the local energy. The angle of rotation can be chosen randomly for each 
spin. However, in order to have minimal autocorrelation times, a reflection — i.e., a rotation 
of a spin by 180 degrees — turns out to be the most efficient overrelaxation move. In one 
overrelaxation sweep this update is applied in sequence to every spin of the lattice. 

Typically, a hybrid Monte-Carlo step consists of two sweeps of the whole lattice via the 
Metropolis algorithm (M), four sweeps of the whole lattice by means of the overrelaxation 
algorithm (O) described above, and four single cluster updates according to the Wolff algo- 
rithm (C). The individual updates are mixed automatically in the program so as to generate 
the update sequence MOCOCMOCOC.As our random number generator, we have uti- 
lized the shift-register generator R1279 defined by the recursion relation X n = X„_ p © X n ~ q 
for (p, q) = (1279, 1063). Generators like this one are among the fastest available. However, 
they are known to cause systematic errors in combination with the Wolff algorithm.^ For 
lags (p, q) as large as the ones used here, these errors are much smaller than typical statistical 
errors. The hybrid nature of the algorithm reduces them further.il 

Using this Monte-Carlo scheme, we have investigated lattice sizes L between L = 12 
and L — 72. The integrated correlation time of the hybrid algorithm is determined by the 
autocorrelation function of the energy or, equivalently, by the autocorrelation function of 
the modulus of the magnetization. Either quantity is 0(3) symmetric, and for sufficiently 
long times, the decay of the corresponding autocorrelation functions is governed by the 
same autocorrelation time. This time scale characterizes the slowest mode of the Wolff 
algorithm, so it also determines the correlation time of our hybrid Monte-Carlo algorithm. 
Note that the autocorrelation time of the Metropolis algorithm is determined by the decay 
of the autocorrelation function of the order parameter, which decays particularly slowly 
near the critical point (critical slowing down). For the hybrid scheme described above, the 
autocorrelation time does not exceed 10 hybrid Monte-Carlo steps for the largest lattice size 
at T = T c . In order not to obtain too strongly correlated initial conditions for the equation 
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of motion, an initial condition is generated every tenth hybrid Monte-Carlo step. 

The integration procedure for the equations of motion is completely separated from the 
Monte-Carlo part of the program. The second-order sublattice decomposition integrator 
described in Ref. EU] is used here. Long-time stability is provided by the exact conservation 
of energy [see Eq. (^j)] and spin normalization. The magnetization is only conserved within 
the accuracy of the discretization, i.e., to second order in the time step. Typical time 
steps St used here range from St = 0.04/ J to St = 0.08/ J, depending on the size of the 
system. For the largest system (L = 72), the total integration time is 77 = 8192/ J; in 
this case St = 0.08/ J was used. Note that the decomposition integrator I (St) has the 
exact time inversion property I (—St) = I^ 1 (St). This guarantees that the time evolution of 
discretization errors, such as those affecting the conservation of the magnetization, does not 
contain systematic drifts.il 

If the algorithm is implemented on a single processor machine, the major part of the CPU 
time is consumed by the integration of the equations of motion. This fraction increases with 
system size because the CPU time needed for the integration scales as 77L 3 , whereas that of 
a hybrid Monte-Carlo step scales as L 3 . If Wolff updates are used exclusively, the average 
scaling is reduced to L 2-r? .il For the purposes of the present investigation the integration 
time 77 has to be chosen such that the slowest spin wave or spin diffusion modes in the system 
can be identified. At T = T c this means that 77 ~ L 3 [see Eq. ([Hp]. Below T c , one must 
have 77 ~ L 2 for an isotropic ferromagnet in order to resolve the slowest spin-wave mode. 
Above T c , the dynamics is dominated by spin diffusion, which also requires 77 ~ L 2 for the 
resolution of the slowest modes. It is therefore very desirable to distribute the integration 
task of the simulation over several processors on a parallel machine. 

A simple and very efficient implementation on a parallel machine with at least four 
processors can be constructed according to the following master-slave idea: The master 
process runs the Monte-Carlo part of the simulation to provide initial conditions on demand, 
and the slave processes integrate the equations of motion for different initial conditions in 
parallel and analyze the correlations. The amount of communication among the processes 
is determined by the transfer of initial configurations from the master to the slaves at the 
beginning of the simulation, and by the transfer of the correlation data from the slaves 
to the master for the final evaluation and output. If N processors in parallel are used 
in this way, the speedup is very close to the theoretical limit iV — 1 for sufficiently large 
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integration times tj. We have implemented the master-slave version of the spin dynamics 
algorithm on the Alpha Linux cluster ALiCE at the Institut fur angewandte Informatik at 
the Bergische Universitat Wuppertal, using up to 32 processors in parallel for the largest 
systems. Communication between the processors is facilitated by the MPI (message passing 
interface) communication library. 

A well-known major problem one is faced with in any computer simulation study of 
critical behavior is how to extract the asymptotic critical behavior from the data. This is 
particularly challenging in our case since we have to cope with two additional complications: 
surface critical behavior and dynamics. In the asymptotic critical regime the value of the 
ratio ri = J\/J does not matter if d = 3 because surfaces of three-dimensional isotropic 
Heisenberg ferromagnets are always disordered in the absence of external fields. Thus such 
systems always belong to the ordinary surface universality class. However, to what extent 
the asymptotic scaling can actually be observed in a computer simulation on finite systems 
is a completely different issue. 

The experience made in a previous study of the static case by one of us@ suggests that 
it should be possible to avoid extended crossover regimes by a careful choice of the ratio 
r\. In order to find out which value of r± is optimal in the sense of giving the largest 
asymptotic regime, we proceed as in Ref. [JT]: We consider the magnetization profile in 
thermal equilibrium, determine r\ in such a way that a discrete version of the Dirichlet 
boundary condition holds, and then verify that this choice suppresses corrections to the 
asymptotic behavior, making the asymptotic regime larger than for alternative values of r%. 
Let us explain this in detail. The equilibrium profile is 



where M tot = ^\ Si is the total magnetization, while z = z 3 + 1/2 with i 3 — 0, . . . , L — 1 
indicates a lattice plane parallel to the surfaces [cf. Eq. ([!])]. Note that according to this 
definition of z, the 'boundary planes' z — 1/2 and z = L — 1/2 of the system are located 
half a lattice constant away from the first and last lattice layers, respectively. With this 



the Ginzburg-Landau Hamiltonian @, where the order parameter of the numerical cell i is 
represented by the spin Si at its center. 

The bulk magnetization can be approximated by the value of the magnetization in the 




(116) 



convention, the lattice model defined by Eq. (P may be viewed as a discrete version of 
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center layer of the system, i.e., = m(z m ia) with z mi< i = L/2 when the number of layers, L, 
is odd. For even L, z mi ^i = (L — l)/2 and z m id,2 = (L + l)/2 are equivalent choices for z mic i; 
in this case, m b = (m(z mid l ) + m(z mid 2 ))/2 is interpreted as the bulk magnetization. The 
ratio m{z)/m^ then depends only on z/L, which motivates us to define the scaling function 

M(z/L) = m(z)/m h , (117) 

where T = T c is assumed. The analysis of the data reveals that the scaling function M(( = 
z/L) can be represented by the simple fit formula 

M(C) = B M [(C + Co)(l - C + (o)}^-™" (H8) 

to a remarkable accuracy. Here Co — Zq/L is the scaled extrapolation parameter. In analyzing 
the data we have accepted the estimates (3 = 0.3662 ± 0.0025 and v = 0.7073 ± 0.0035 of 
Ref. EH, and utilized the value ft = 0.834(6) of Ref. EH. 



From a least square fit of Eq. ( |118| ) to the data for various system sizes we obtain the 
extrapolation parameter zq in units of the lattice spacing. For the choices J±/J = 0.3 
and 1.0, we find z ~ —0.34 and z ~ 0.46, respectively. On the other hand, z vanishes 
within the statistical errors if Ji/J = 0.73.0 To put these findings in perspective, some 



explanations are necessary. Owing to our definition of z [given below Eq. (|116 )1, a fit curve 
( |118|) with scaled extrapolation parameter (q = means that the measured magnetization 
profile extrapolates to zero half a lattice constant away from the outermost layers is = 
and «3 = L — 1 of our lattice model. In this sense, the profile satisfies a Dirichlet boundary 
condition on the scale of the lattice constant in this special case. 

Let us emphasize that such a boundary condition on a microscopic scale must not be 
confused with the Dirichlet boundary condition which the order parameter satisfies at the 
ordinary transition on long scales, irrespective of the precise value of the ratio r\ = J\j J ■ 
The latter is an asymptotic long-scale "property, associated with the corresponding ordinary 
fixed point of the RG, and hence universal. By contrast, the boundary condition that the 
order-parameter profile of a given microscopic model is found to satisfy on microscopic scales 
generally depends on microscopic details, and is therefore a nonuniversal property (cf. the 
discussion in Sec. III.C.9 of Ref. ||). 

On the level of a mesoscopic description via our continuum model, a Dirichlet boundary 
condition can be enforced on the mesoscopic length scales on which such a description 
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makes sense (several lattice constants) by setting the enhancement variable Co to the fixed- 
point value +00. For values cq < 00, the Dirichlet boundary condition does not hold on 
mesoscopic scales, neither for the regularized nor for the renormalized theory. In other words, 
a Co-dependent extrapolation parameter z 7^ occurs. This deviation from the Dirichlet 
boundary condition corresponds to a correction to scaling: It is irrelevant inasmuch as it 
vanishes in the limit zq/z — > 0. Choosing a particular value r± = rp for the ratio of 
interaction constants such that zq ~ is an appealing way to mimic the Dirichlet boundary 
condition of the Cq = 00 continuum theory on a lattice. As we know already for the static 



case from Ref. [31|, and will verify for the dynamic theory below, this choice of r\ suppresses 
corrections to scaling and hence enlarges the regime in which the asymptotic scaling behavior 
is observed. 

The optimal value which yields a vanishing extrapolation parameter zq for temperatures 
sufficiently close to T c and moderately large lattice dimensions L is r[ D ^ = 0.73. For values 
close to this optimal one, we have 

z = a (n - r[ D) ) + 0[(n - r[ D) ) 2 } , (119) 

where a is a factor of order unity. This behavior was already obtained in Ref. |HT|, where 
crossover effects and the behavior of the order parameter profile as a function of r\ were 
investigated in more detail for the static case. The present work confirms these findings: 
Our results for the dynamic surface structure function presented in the next section are fully 
consistent with them. 



V. THE DYNAMIC SURFACE STRUCTURE FUNCTION 

Our simulation results are displayed in Figs. They were obtained for T = T c , 

L = 72, and the total integration time 77 = 8192/ J. For the smallest accessible frequency 
^min = 2tt/ti, finite-size effects turned out to be negligible. 

In Fig. [I] the structure function Cii(0, u) is shown for different values of r\ in comparison 
with Eq. (|115|) . The exponent (3 + 1 — ?7j| rd )/3 has the value 0.856 ± 0.005 that follows from 
the estimate a(d=3) = 2.482 ± 0.002 obtained by the substitution i](d=3) = 0.036 ± 0.0040 
in Eq. (|7T|) and the current estimate of the surface correlation exponent ?7j| rd (<i=3) = 1.358 ± 
0.0120 of the ordinary transition. 
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FIG. 1: Structure function C n (0,u) for n = J x /J = 0.3 (x), 0.73 (+), and 1.0 (*). Error bars (one 
standard deviation) are smaller than the symbol sizes. The solid lines indicate the theoretically 
expected power law (|115| ) for uj — > 0. 



The dependence of Cu(0,u) on n is particularly interesting. If n is small (r x = 0.3, x), 
our simulation data approach the asymptotic power law ( |115| ) from above, whereas for larger 
values of r% (ri = 1, *), the asymptotic power law is approached from below. In the latter 
case, the data even remain outside the asymptotic regime for the frequency range shown 
in Fig. |l|. The best agreement with Eq. ( |115| ) over the largest frequency range is obtained 
for the choice T\ = 0.73(+), which has already been identified as optimal in the sense that 
the extrapolation parameter z for the magnetization profile vanishes (see Eq. ( |119| )). The 
deviations from the power law ( 115 ) for r\ ^ 0.73 can apparently be attributed to dynamic 
surface-induced corrections to the asymptotic behavior that originate from the nonzero, 
rx-dependent value of the extrapolation parameter z . 

Fig. H shows a scaling plot of Cn(p, uj), where p = (^|, 0, 0) is oriented along the surface. 
As one sees, the scaling regime in x shrinks as the mode index n is increased from l(x) to 
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FIG. 2: Scaling function o{x) according to ( 114 ). Data obtained for r\ = J\/J = 0.73 and 
p = (§f , 0, 0), with n = 1, . . . ,5, are shown. Error bars (one standard deviation) are smaller than 
the symbol sizes. The data follow Eq. (|120| ) up to x ~ 1. The data for x > 1 are outside the scaling 
regime. 



5(b). For x < 1, the shape of the scaling function in Eq. (|114|) is described surprisingly well 
by the fit functionlli 



a (x) 



a [l + (x/x o r]^ d -^ , ( 120 ) 

which is inspired by the known zero- loop result Bi'i The exponent at the square bracket is 
chosen so as to reproduce Eq. ( |115|) in the limit x — > oo (p — > at fixed uj ^ 0). The 
amplitude <7q and the crossover parameter xq are used as fit parameters. 

The agreement between the data displayed in Figs. [I] and ||| and the scaling forms (|114|) 
and ( 115 ) is quite satisfactory. However, on closer inspection small deviations are found to 
remain. Note that as pointed out at the end of the previous section and in analogy with 



the results of Ref. [31] for the equilibrium case, the choice r\ = 0.73 yields indeed the largest 
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regime in which asymptotic scaling holds (see Fig. [I]). 
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FIG. 3: Scaling plot of the surface structure function for n = J\/ J = 0.3 and p = (^,0,0), with 
n = 1, . . . ,5. Error bars (one standard deviation) are smaller than the symbol sizes. The data do 
not follow Eq. (|T20|). 



Hence we expect the choice r% = 0.73 to be optimal also for the surface structure factor 
at finite momentum transfer p. Our results for r\ = 0.3 and r± — 1.0 depicted in Figs. ||] and 
£| confirm this expectation. Figure ||] shows that the data for r\ = 0.3 approach the scaling 
function a(x) of Fig. ^| monotonically from below as the mode number n is increased from 1 
to 5. The nonasymptotic surface- induced corrections are so large that the data for different 
n (i.e., momentum transfers p) are well separated even on a logarithmic scale. In other 
words, no data collapse nearly as nice as in Fig. |2| occurs, although the crossover parameter 
x appears to be consistent with the results displayed there. 

Our results for r\ = 1.0 (see Fig. ^) show a similar behavior, except that the scaling 
function cr(x) now is monotonically approached from above as the mode number n increases. 
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The crossover parameter xq is again consistent with our findings in Figs. and [| The 
nonasymptotic surface-induced corrections are as large as in Fig. [3] but have opposite sign. 
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FIG. 4: Scaling plot of the surface structure function for r\ = 1.0 and p = (^,0,0), with n = 
1, . . . ,5. Error bars (one standard deviation) are smaller than the symbol sizes. The data do not 
follow Eq. (pop. 



From these results a consistent numerical picture of dynamic surface scaling emerges. For 
the optimal choice r% = 0.73, our simulation data for the dynamic structure function bear 
out quite convincingly the asymptotic behavior we predicted on the basis of our RG work. 
For values of r\ deviating significantly from 0.73, the data exhibit pronounced corrections 
to scaling. These entail that our data for values like r\ — 1 and 0.3 do not exhibit directly 
the asymptotic scaling form and power law of the surface structure function Cn(p,uj) for 
nonzero and zero parallel momentum p, respectively. Yet they seem to be consistent both 
with the theoretically predicted asymptotic behavior as well as with the one extracted from 
our simulation data for n = 0.73 because the observed deviations appear to be attributable 
to corrections to scaling. 
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In order to demonstrate universality, it would clearly be very desirable to verify the 
approach to the asymptotic critical behavior also for unfavorable values of r\ like 1.0 and 
0.3. One conceivable way of trying to reach this goal is by means of brute-force numerical 
means. However, in view of the enormous numerical effort that was necessary to produce 
the above simulation results, we do not consider this to be promising strategy at present. 

We believe that a better strategy would be the incorporation of nonasymptotic correction 
terms in the analysis of the simulation data. Unfortunately, there are various sources of such 
corrections, and detailed knowledge about their form and strength is either scarce or not 
available. A systematic investigation of the various kinds of nonasymptotic corrections of 
static and dynamic origin that might play a role in the analysis of the surface critical behavior 
we are concerned with here evidently requires more numerical and analytical work, and is 
beyond the scope of this paper. Let us therefore simply discuss some possible sources of 
deviations from scaling, beginning with the ones that do not correspond to corrections to 
scaling. 

Two obvious sources of this latter kind are insufficient momentum and frequency resolu- 
tion. By virtue of the relation Sp = 2tc / L the momentum resolution Sp is intimately linked 
to the system size L, which despite formidable progress in simulation techniques and com- 
putational resources still is a serious limiting factor. The frequency resolution Su = 2it / tj is 
limited by the total integration time 77. From our data for Cu(p, t) (not shown) we conclude 
that 77 is sufficiently long. The frequency resolution Suj/J ~ 7.7 x 10~ 4 that is available 
here rivals that of neutron scattering experiments.0 The momentum resolution is given by 
Sp ~ 0.087 in units of the inverse lattice constant for our largest systems, and is therefore 
much more restrictive. 

One familiar type of corrections to scaling are those induced by deviations of the coupling 
constant u from its fixed point value u* . Just as in the static case, they are governed by 
the Wegner exponent u u = j3' u (u*) whose value is ~ 0.8 in d = 3 dimensions. Analogous 
corrections to scaling result from the RG flow of the mode-coupling interaction constant 
/. Upon linearizing the flow about the infrared-stable fixed point (u*,f*), one obtains in 
addition to uj u a second correction-to-scaling exponent, uif = f3'Ju*,f*), which in contrast 
to the former is of purely dynamic origin. We are not aware of any reliable estimates of Uf 
for d = 3.i 

Other potentially dangerous corrections might be caused by a previously mentioned im- 
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portant difference of the dynamics of the simulated lattice model and the semi-infinite con- 
tinuum model J: that the former conserves the energy while the latter does not. In Appendix 
[A| we generalize model J for the bulk case by incorporating energy conservation. The result- 
ing model is analogous to model D (with n = 3 components), and reduces to this when the 
mode-coupling interaction constant / is set to zero. We show that the ratio A/A^- 1 of trans- 
port coefficients (where \^ denotes the analog of A for the energy density 8) transforms 
under the RG as £ 3 ~ 2 in the long length-scale limit i — > 0. Since 3 — 2 = (d — 2 + i])/2, which 
is ~ 0.5 in three dimensions, the ratio approaches indeed zero, albeit with a considerably 
smaller power than in the case of model D (where the value of this exponent is 2 — 77 ~ 2). 

Hence two conclusions may be drawn: First, in order to obtain the asymptotic critical 
behavior of order parameter cumulants we can take the limit A/A (£) -> 0. If one sets 1/A< f > = 
from the outset, the energy density relaxes instantaneously, This means that the effects 
produced by the coupling to the energy density correspond to a change of the parameters 
of the original model J, up to corrections due to irrelevant operators. In other words, the 
energy conservation should not affect the asymptotic critical behavior, so that our lattice 
model should belong to the universality class of our semi-infinite model J. Second, we cannot 
rule out that the corrections to the asymptotic behavior induced by the conservation of the 
energy are less important for an improved analysis of the numerical data presented above 
than the previously mentioned corrections to scaling. To assess the relative importance of the 
various types of corrections to scaling seems difficult without reliable additional information 
based on detailed calculations. 

The corrections to scaling we have just considered are associated with irrelevant bulk 
variables and hence are not specific to systems with surfaces. Analogous ones are induced 
by irrelevant surface variables. A well-known example are the corrections ~ I resulting from 
deviations of 1/cq (the reciprocal surface enhancement variable) from the fixed-point value 
l/c* rd = 0. One evident consequence of such corrections (which is, however, not the only 
one when Landau theory is not valid) is that the Dirichlet boundary condition the order- 
parameter density satisfies at the ordinary fixed point ceases to hold. In view of the two 
observations made above — namely, (i) that the choice T\ = 0.73 suppresses corrections to 
scaling both in the case of the dynamic structure functions as well as in static quantities,^ 
and (ii) that the deviations from scaling according to Figs. ^] and ^ have opposite signs 
depending on whether 7*1 is bigger or smaller than the optimal value 0.73 — the corrections 
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to scaling which the finiteness of Co induces appear to play a major role. 

VI. SUMMARY AND CONCLUSIONS 

We have presented a detailed study of the surface critical behavior of isotropic Heisenberg 
ferromagnets at the ordinary transition, using both sophisticated analytical tools as well 
as high-precision Monte-Carlo spin dynamics simulations. In the analytical part of our 
work a continuum model that represents the corresponding universality class of bulk and 
surface critical behavior — namely, an appropriate semi-infinite extension of model J — has 
been formulated and its field theory constructed. To this end we have determined the 
relevant boundary contributions to the dynamic action functional which are compatible 
with the general features (symmetries, detailed balance, locality assumptions, conservation 
laws, etc.) of this class of systems. These were shown to correspond to boundary conditions 
for the resulting dynamic field theory. 

In order to investigate the critical behavior of this semi-infinite model J, we have employed 
two distinct RG schemes. The first is a massless one based on the expansion about the 
upper critical dimension dj — 6. To avoid the difficulties it has in handling the problem 
adequately below the upper critical dimension d* = 4 of the static theory, we have designed 
and utilized an appropriate extension of the massive RG scheme for bounded systems of 
Diehl and Shpot. As usual, this works in fixed dimensions, and avoids the dimensionality 
expansion. By combining the resulting RG equations with the boundary operator expansion 
(ff5|), we have been able to obtain detailed predictions for the scaling behavior of the surface 
structure function (|5]). The involved critical exponents which govern power laws like (|115| ) 
are related to known static bulk and surface critical exponents. In particular, there is no 
independent new dynamic exponent associated with the surface.il 

We have checked our predictions by means of extensive Monte-Carlo spin dynamics sim- 
ulations. Our results depicted in Figs. [I] and Q corroborate the predicted dynamic scaling 
behavior. In order to obtain such a good manifestation of scaling, we have found it helpful 
and necessary to choose the ratio r% = Ji/Joi the surface and bulk exchange integrals J\ 
and J such that corrections to scaling are suppressed. To achieve this objective we have 
optimized the value of r\ by requiring that the Monte Carlo data yield an equilibrium order- 
parameter profile which satisfies a Dirichlet boundary condition on the scale of the lattice 
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constant in the sense described in Sec. |V|. 

According to our results displayed in Figs. |I|-f|, this procedure is quite effective: For 
ri = 0.73 and values inside a narrow range around this optimal one, the simulation data for 
the dynamic structure functions exhibit clear evidence of dynamic scaling in conformity with 
our predictions. On the other hand, for values of r\ outside this regime, the data collapse 
is poor and the asymptotic behavior cannot be identified in a convincing fashion. These 
observations are in conformity with those made in a previous Monte Carlo investigation of 
the static surface critical behavior by one of us.0 However, in the study of static quantities 
one is in a much better position because the scaling regimes can be reached reasonably well 
even for non-optimal values of r\. This fact lends support to our belief that the dynamic 
scaling behavior seen for r\ = 0.73 can be trusted. 

Finally, let us note that we have not taken into account any dipolar forces in our study. 
To check our result by experiments (as should become feasible in the near future thanks to 
facilities like the X-ray free electron lasei0), one must choose systems for which such forces 
are negligible. Since even weak dipolar forces lead to the formation of domains, one must 
also make sure that single domains are investigated. 
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APPENDIX A: THE IRRELEVANCE OF ENERGY CONSERVATION 

The dynamics of the lattice model we simulated, defined via the equations of motion 
(@), conserves the energy. Our aim here is to show that this feature does not change the 
universality class, which should therefore be represented by the semi-infinite model J em- 
ployed in our RG study. For the sake of simplicity, we restrict ourselves to demonstrating 
the irrelevance of energy conservation for the bulk case. The extension to the semi-infinite 
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case should be obvious. 

Conservation of energy means that the energy density £(x, t) is a conserved density and 
hence a slow variable which should be retained in a coarse-grained description on mesoscopic 
time scales. Now, for vanishing mode-coupling constant / , model J reduces to model B. 
How to modify the latter to account for energy conservation is well known and leads to 
model D.i We can adapt the dynamics of model J along similar lines. The obvious result 
is a modification of model J that differs from model D through the addition of the former's 
mode-coupling terms. The Langevin equations of this two-density model, which we call J', 
read 

t) = A (a ^ + f ^ x 4^ + C(x, t) (Al) 



and 



S(x,t) = \^A 5 -^ + #(x,t), (A2) 



where 



i ( v0) 2 + ^ 2 + ^i0r + ^ 2 +|^ 2 



(A3) 



is the familiar Hamiltonian employed in the definition of models C, D, and E (for the here 
considered case of an n = 3 component order parameter <fi) . Both £ as well as d are Gaussian 
random forces with zero average; their variances are given by Eq. (|8]) and 

{&(x, t) 0(x', f)) = -2X^ ] A 5{x-x') 8{t-t') , (A4) 

respectively. 

In the absence of coupling between the order parameter <f> and the energy density £ , i.e., 
for 7q = 0, the dynamic exponent of £ takes its Gaussian value 

U = 2 , (A5) 

corresponding to ordinary diffusion. It is not difficult to see that this result remains valid 
for 7q 7^ 0. For d > 4, this follows immediately from the observation that 70 is irrelevant in 
the RG sense. 

In studying the more interesting case d < 4, we can benefit from well-known results 
for the static theory described by the Hamiltonian (A"3), which is equivalent to the 4 - 



Hamiltonian (|6|) except for a change u — » U = u — 3 7q of the interaction constant (see, 
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for example, Refs. [59| and Ref. ^ and their references). As a consequence, the corresponding 
static renormalization functions can be expressed in terms of those of the </> 4 -theory. In 
particular, the renormalized analog U of U can be introduced in analogy to Eq. ( |79|) via 
Uq = Z U (U) m A ~ d U , where Z u (u) is the renormalization function of Sec. [111D 3| . Likewise, 
the renormalization factor Z 1 (U, 7 2 ) which relates to its renormalized 

counterpart 7 can be expressed in terms of known renormalization functions of the 4 - 
theory. (It is a product of Z^iJJ) and a factor linear in 7 2 whose ^-dependent coefficient 
derives from the additive counterterm that the 4 -vertex function Ta^a? requires. 

The resulting RG flow of U and 7 has two nontrivial fixed points at U — u*: one at 
7* = 0, and another one at (7*) 2 = const a/z/. The slopes d(3r y 2{U=u* , ■y 2 )/d r y 2 of the beta 
function /3 7 2 = rnd m \ -y 2 at these fixed points are given by —ajv and ajv, respectively. 

Since a < in the three-component case we are concerned with (a ~ —0.12 for d = 3, 
according to Ref. |56|), the infrared-stable fixed point is (U, 7) = {u* , 0). The results of Ref. ^9 
imply that the running coupling constant associated with 7 tends to zero as m~ a ^ 2u in the 
limit m — > 0. Thus the energy density decouples asymptotically from the order parameter, 
so that the result (|A5|) applies. 

We can introduce the renormalized transport coefficient via = vrC 2 Zg\^ £ \ 
where Zg, the static renormalization factor of the energy density, takes the value 1 at the 
infrared-stable fixed point. The ratio of transport coefficients A/A^ has the asymptotic 
scale dependence ~ m?~ u . If we substitute the values ([71]) and ( |A5| ) for 3 and 3^, the 
exponent becomes 3— is — (c? — 2 — ?7)/2. Since this is positive in three dimensions, the ratio 
approaches zero in the long-scale limit m — > 0. The upshot is that the critical dynamics of 
the order parameter remains unaffected by the coupling to the energy density, as claimed. 



* Present address: Institut fur Theoretische und Angewandte Physik, Universitat Stuttgart, 70550 
Stuttgart and Max-Planck-Institut fiir Metallforschung, Heisenbergstr. 1, 70569 Stuttgart, Fed- 
eral Republic of Germany 
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